home *** CD-ROM | disk | FTP | other *** search
/ Computer Select (Limited Edition) / Computer Select.iso / dobbs / v17n04 / wavelet.asc < prev    next >
Encoding:
Text File  |  1992-03-10  |  11.0 KB  |  328 lines
  1. _THE FAST WAVELET TRANSFORM_
  2. by Mac A. Cody
  3.  
  4. [LISTING ONE]
  5.  
  6. #define WAVE_MGT
  7. #include <alloc.h>
  8. #include "wave_mgt.h"
  9. double **BuildTreeStorage(int inlength, int levels)
  10. {
  11.   double **tree;
  12.   int i, j;
  13.   /* create decomposition tree */
  14.   tree = (double **) calloc(2 * levels, sizeof(double *));
  15.   j = inlength;
  16.   for (i = 0; i < levels; i++)
  17.   {
  18.     j /= 2;
  19.     if (j == 0)
  20.     {
  21.       levels = i;
  22. /*    printf("\nToo many levels requested for available data\n");
  23.       printf("Number of transform levels now set to %d\n", levels); */
  24.       continue;
  25.     }
  26.     tree[2 * i] = (double *) calloc((j + 5), sizeof(double));
  27.     tree[2 * i + 1] = (double *) calloc((j + 5), sizeof(double));
  28.   }
  29.   return tree;
  30. }
  31. void DestroyTreeStorage(double **tree, int levels)
  32. {
  33.   char i;
  34.   for (i = (2 * levels - 1); i >= 0; i--)
  35.     free(tree[i]);
  36.   free(tree);
  37. }
  38. void TreeCopy(double **TreeDest, double **TreeSrc, int siglen, int levels)
  39. {
  40.   int i, j;
  41.   for (i = 0; i < levels; i++)
  42.   {
  43.     siglen /= 2;
  44.     for (j = 0; j < siglen + 5; j++)
  45.     {
  46.       if ((i + 1) == levels)
  47.         TreeDest[2 * i][j] = TreeSrc[2 * i][j];
  48.       else
  49.         TreeDest[2 * i][j] = 0.0;
  50.       TreeDest[(2 * i) + 1][j] = TreeSrc[(2 * i) + 1][j];
  51.     }
  52.   }
  53. }
  54. void TreeZero(double **Tree, int siglen, int levels)
  55. {
  56.   int i, j;
  57.   for (i = 0; i < levels; i++)
  58.   {
  59.     siglen /= 2;
  60.     for (j = 0; j < siglen + 5; j++)
  61.     {
  62.       Tree[2 * i][j] = 0.0;
  63.       Tree[(2 * i) + 1][j] = 0.0;
  64.     }
  65.   }
  66. }
  67. void ZeroTreeDetail(double **Tree, int siglen, int levels)
  68. {
  69.   int i, j;
  70.   for (i = 0; i < levels; i++)
  71.   {
  72.     siglen /= 2;
  73.     for (j = 0; j < siglen + 5; j++)
  74.       Tree[(2 * i) + 1][j] = 0.0;
  75.   }
  76. }
  77.  
  78.  
  79. [LISTING TWO]
  80.  
  81. /* WAVELET.C */
  82. #include <math.h>
  83. typedef enum { DECOMP, RECON } wavetype;
  84. #include "wavelet.h"
  85. void WaveletCoeffs(double alpha, double beta, double *wavecoeffs)
  86. {
  87.   double tcosa, tcosb, tsina, tsinb;
  88.   char i;
  89.   /* precalculate cosine of alpha and sine of beta to reduce */
  90.   /* processing time */
  91.   tcosa = cos(alpha);
  92.   tcosb = cos(beta);
  93.   tsina = sin(alpha);
  94.   tsinb = sin(beta);
  95.   /* calculate first two wavelet coefficients,  a = a(-2) and b = a(-1) */
  96.   wavecoeffs[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
  97.                                       + 2.0 * tsinb * tcosa) / 4.0;
  98.   wavecoeffs[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
  99.                                       - 2.0 * tsinb * tcosa) / 4.0;
  100.   /* precalculate cosine and sine of alpha minus beta to reduce */
  101.   /* processing time */
  102.   tcosa = cos(alpha - beta);
  103.   tsina = sin(alpha - beta);
  104.   /* calculate last four wavelet coefficients  c = a(0), d = a(1), */
  105.   /* e = a(2), and f = a(3) */
  106.   wavecoeffs[2]  = (1.0 + tcosa + tsina) / 2.0;
  107.   wavecoeffs[3]  = (1.0 + tcosa - tsina) / 2.0;
  108.   wavecoeffs[4]  = 1 - wavecoeffs[0] - wavecoeffs[2];
  109.   wavecoeffs[5]  = 1 - wavecoeffs[1] - wavecoeffs[3];
  110.   /* zero out very small coefficient values caused by truncation error */
  111.   for (i = 0; i < 6; i++)
  112.   {
  113.     if (fabs(wavecoeffs[i]) < 1.0e-15)
  114.       wavecoeffs[i] = 0.0;
  115.   }
  116. }
  117. char MakeWaveletFilters(double *wavecoeffs, double *Lfilter,
  118.                                double *Hfilter, wavetype transform)
  119. {
  120.   char i, j, k, filterlength;
  121.   /* find the first non-zero wavelet coefficient */
  122.   i = 0;
  123.   while(wavecoeffs[i] == 0.0)
  124.     i++;
  125.   /* find the last non-zero wavelet coefficient */
  126.   j = 5;
  127.   while(wavecoeffs[j] == 0.0)
  128.     j--;
  129.  /* form decomposition filters h~ and g~ or reconstruction filters h and g.
  130.  Division by 2 in construction of decomposition filters is for normalization */
  131.   filterlength = j - i + 1;
  132.   for(k = 0; k < filterlength; k++)
  133.   {
  134.     if (transform == DECOMP)
  135.     {
  136.       Lfilter[k] = wavecoeffs[j--] / 2.0;
  137.       Hfilter[k] = (double) (((i++ & 0x01) * 2) - 1) * wavecoeffs[i] / 2.0;
  138.     }
  139.     else
  140.     {
  141.       Lfilter[k] = wavecoeffs[i++];
  142.       Hfilter[k] = (double) (((j-- & 0x01) * 2) - 1) * wavecoeffs[j];
  143.     }
  144.   }
  145.   /* clear out the additional array locations, if any */
  146.   while (k < 6)
  147.   {
  148.     Lfilter[k] = 0.0;
  149.     Hfilter[k++] = 0.0;
  150.   }
  151.   return filterlength;
  152. }
  153. double DotP(double *data, double *filter, char filtlen)
  154. {
  155.   char i;
  156.   double sum;
  157.   sum = 0.0;
  158.   for (i = 0; i < filtlen; i++)
  159.     sum += *data-- * *filter++; /* decreasing data pointer is */
  160.                                 /* moving backward in time */
  161.   return sum;
  162. }
  163. void ConvolveDec2(double *input_sequence, int inp_length,
  164.           double *filter, char filtlen, double *output_sequence)
  165. /* convolve the input sequence with the filter and decimate by two */
  166. {
  167.  int i;
  168.  char shortlen, offset;
  169.  for(i = 0; (i <= inp_length + 8) && ((i - filtlen) <= inp_length + 8); i += 2)
  170.   {
  171.     if (i < filtlen)
  172.       *output_sequence++ = DotP(input_sequence + i, filter, i + 1);
  173.     else if (i > (inp_length + 5))
  174.     {
  175.       shortlen = filtlen - (char) (i - inp_length - 4);
  176.       offset = (char) (i - inp_length - 4);
  177.       *output_sequence++ = DotP(input_sequence + inp_length + 4,
  178.                                         filter + offset, shortlen); 
  179.     }
  180.     else
  181.       *output_sequence++ = DotP(input_sequence + i, filter, filtlen);
  182.   }
  183. }
  184. int DecomposeBranches(double *In, int Inlen, double *Lfilter,
  185.     double *Hfilter, char filtlen, double *OutL, double *OutH)
  186.   /* Take input data and filters and form two branches of have the
  187.      original length. Length of branches is returned */
  188. {
  189.   ConvolveDec2(In, Inlen, Lfilter, filtlen, OutL);
  190.   ConvolveDec2(In, Inlen, Hfilter, filtlen, OutH);
  191.   return (Inlen / 2);
  192. }
  193. void WaveletDecomposition(double *InData, int Inlength, double *Lfilter,
  194.       double *Hfilter, char filtlen, char levels, double **OutData)
  195. /* Assumes input data has 2 ^ (levels) data points/unit interval. First InData 
  196. is input signal; all others are intermediate approximation coefficients */
  197. {
  198.   char i;
  199.   for (i = 0; i < levels; i++)
  200.   {
  201.     Inlength = DecomposeBranches(InData, Inlength, Lfilter, Hfilter,
  202.                         filtlen, OutData[2 * i], OutData[(2 * i) + 1]);
  203.     InData = OutData[2 * i];
  204.   }
  205. }
  206. double DotpEven(double *data, double *filter, char filtlen)
  207. {
  208.   char i;
  209.   double sum;
  210.   sum = 0.0;
  211.   for (i = 0; i < filtlen; i += 2)
  212.     sum += *data-- * filter[i]; /* decreasing data pointer is moving */
  213.                                 /* backward in time */
  214.   return sum;
  215. }
  216. double DotpOdd(double *data, double *filter, char filtlen)
  217. {
  218.   char i;
  219.   double sum;
  220.   sum = 0.0;
  221.   for (i = 1; i < filtlen; i += 2)
  222.     sum += *data-- * filter[i];  /* decreasing data pointer is moving */
  223.                                  /* backward in time */
  224.   return sum;
  225. }
  226. void ConvolveInt2(double *input_sequence, int inp_length, double *filter,
  227.           char filtlen, char sum_output, double *output_sequence)
  228.   /* insert zeros between each element of the input sequence and
  229.      convolve with the filter to interpolate the data */
  230. {
  231.   int i;
  232.   if (sum_output) /* summation with previous convolution if true */
  233.   {
  234.       /* every other dot product interpolates the data */
  235.     for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
  236.     {
  237.       *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
  238.       *output_sequence++ += DotpEven(input_sequence + i + 1, filter, filtlen);
  239.     }
  240.     *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
  241.   }
  242.   else /* first convolution of pair if false */
  243.   {
  244.       /* every other dot product interpolates the data */
  245.     for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
  246.     {
  247.       *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
  248.       *output_sequence++ = DotpEven(input_sequence + i + 1, filter, filtlen);
  249.     }
  250.     *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
  251.   }
  252. }
  253. int ReconstructBranches(double *InL, double *InH, int Inlen,
  254.         double *Lfilter, double *Hfilter, char filtlen, double *Out)
  255.   /* Take input data and filters and form two branches of have 
  256.      original length. length of branches is returned */
  257. {
  258.   ConvolveInt2(InL, Inlen, Lfilter, filtlen, 0, Out);
  259.   ConvolveInt2(InH, Inlen, Hfilter, filtlen, 1, Out);
  260.   return Inlen * 2;
  261. }
  262. void WaveletReconstruction(double **InData, int Inlength, double *Lfilter,
  263.               double *Hfilter, char filtlen, char levels, double *OutData)
  264.   /* assumes that input data has 2 ^ (levels) data points per unit interval */
  265. {
  266.   double *Output;
  267.   char i;
  268.   Inlength = Inlength >> levels;
  269.   /* Destination of all but last branch reconstruction is the next
  270.      higher intermediate approximation */
  271.   for (i = levels - 1; i > 0; i--)
  272.   {
  273.     Output = InData[2 * (i - 1)];
  274.     Inlength = ReconstructBranches(InData[2 * i], InData[(2 * i) + 1],
  275.                           Inlength, Lfilter, Hfilter, filtlen, Output);
  276.   }
  277.   /* Destination of the last branch reconstruction is the output data */
  278.   ReconstructBranches(InData[0], InData[1], Inlength, Lfilter, Hfilter,
  279.                                                       filtlen, OutData);
  280. }
  281. double CalculateMSE(double *DataSet1, double *DataSet2, int length)
  282. {
  283.   /* calculate mean squared error of output of reconstruction as
  284.      compared to the original input data */
  285.   int i;
  286.   double pointdiff, topsum, botsum;
  287.   topsum = botsum = 0.0;
  288.   for (i = 0; i < length; i++)
  289.   {
  290.     pointdiff = DataSet1[i] - DataSet2[i];
  291.     topsum += pointdiff * pointdiff;
  292.     botsum += DataSet1[i] * DataSet1[i];
  293.   }
  294.   return topsum / botsum;
  295. }
  296.  
  297.  
  298. [LISTING THREE]
  299.  
  300. /* WAVE_MGT.H */
  301. double **BuildTreeStorage(int inlength, int levels);
  302. void DestroyTreeStorage(double **tree, int levels);
  303. void TreeCopy(double **TreeDest, double **TreeSrc, int siglen, int levels);
  304. void TreeZero(double **Tree, int siglen, int levels);
  305. void ZeroTreeDetail(double **Tree, int siglen, int levels);
  306. /* WAVELET.H */
  307. void WaveletCoeffs(double alpha, double beta, double *wavecoeffs);
  308. char MakeWaveletFilters(double *wavecoeffs, double *Lfilter,
  309.                 double *Hfilter, wavetype transform);
  310. double DotP(double *data, double *filter, char filtlength);
  311. void ConvolveDec2(double *input_sequence, int inp_length,
  312.         double *filter, char filtlen, double *output_sequence);
  313. int DecomposeBranches(double *In, int Inlen, double *Lfilter,
  314.       double *Hfilter, char filtlen, double *OutL, double *OutH);
  315. void WaveletDecomposition(double *InData, int Inlength, double *Lfilter,
  316.       double *Hfilter, char filtlen, char levels, double **OutData);
  317. double DotpEven(double *data, double *filter, char filtlength);
  318. double DotpOdd(double *data, double *filter, char filtlength);
  319. void ConvolveInt2(double *input_sequence, int inp_length, double *filter,
  320.         char filtlen, char sum_output, double *output_sequence);
  321. int ReconstructBranches(double *InL, double *InH, int Inlen,
  322.       double *Lfilter, double *Hfilter, char filtlen, double *Out);
  323. void WaveletReconstruction(double **InData, int Inlength, double *Lfilter,
  324.       double *Hfilter, char filtlen, char levels, double *OutData);
  325. double CalculateMSE(double *DataSet1, double *DataSet2, int length);
  326.  
  327.  
  328.